[NCLASSIFIED _ 

JR?jY~tLASSIFICATION  OF  THIS  PAGE 


JMENTATION  PAGE 


AD-A211  034 


DECLASSIFICATION  /  DOWNGRADING  SCHEDULE 


Form  Approved 
OMB  No.  0704-0188 


1b  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION /AVAILABILITY  OF 

DISTRIBUTION  UNLIMITED 


Jfi£.  4 


PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

FJSRL-JR-88-0008 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S)  . 

% 

vt 


{1 


NAME  OF  PERFORMING  ORGANIZATION 
Frank  J.  Seiler  Research  Lab 


6b.  OFFICE  SYMBOL 
(If  applicable) 

FJSRL/NC 


7a.  NAME  OF  MONITORING  ORGANIZA; 


ADDRESS  (City,  State,  and  ZIP  Code) 
FJSRL/NC 
USAF  Academy 
Colorado  80840-6528 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 


aTE 


c/rrto 


.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION  AF  Office  of 

Scientific  Research 


8b.  OFFICE  SYMBOL 
(If  applicable) 

AFOSR 


9  PROCUREMENT  INSTRUMENT  IDENTIFI 


N  NUMBER 


Bldg  410 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

DC  20332 

ELEMENT  NO 

NO. 

NO 

ACCESSION  NO 

6U02F 

2303 

F3 

07 

title  (include  Security  Classification)  Numerical  Sensitivity  of  Trajectories  across  Conformational 

Energy  Hypersur faces  from  Geometry  Optimized  Molecular  Orbital  Calculations:  AMI,  JpjJggQ^jjd 


l.  PERSONAL  AUTHOR(S) 

Donald  B.  Boyd,  David  W.  Smith,  James  J.P.  Stewart,  and  Erich  Wimmer 


3a.  TYPE  OF  REPORT 

13b  TIME  COVERED 

14.  DATE  OF  REPORT  (Year.  Month.  Day) 

15.  PAGE  COUNT 

Journal  Article 

FROM  TO 

87/09/22 

12 

S.  SUPPLEMENTARY  NOTATION 


7.  COSATI  CODES  j 

FIELD 

GROUP 

SUB-GROUP 

18  SUBJECT  TERMS  (Continue  on  reverse 


if  necessary  and  identify  by  block  number) 


9.  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

The  monocyclic  B-lactam  ((4(S)-methyl-2-oxo-l-azetidinyl)thia)acetic  acid  was  studied  by  the 
semi-empirical  molecular  orbital  methods  AMI,  MND0,  and  MINDO/3.  Using  the  reaction 
coordinate  option  in  the  program  MOPAC  on  VAX  and  Cray  X-MP  computers,  the  potential  energy 
curve  was  calculated  for  rotation  of  the  C2-NI-S-C  torsional  angle  in  the  conformationally 
flexible  side  chain  while  optimizing  all  other  geometrical  variables  in  the  molecule.  The 
trajectory  taken  during  geometry  optimization  was  found  to  be  sensitive  to  the  computer,  the 
program  version,  the  convergence  criteria,  and  the  degree  of  code  optimization  used  in  the 
calculation.  In  order  to  reduce  the  likelihood  of  spurious  results,  conformational  or 
reaction  energy  hypersurfaces  need  to  be  calculated  with  the  more  precise  SCF  convergence 
and  minimization  criteria  available  in  programs  for  MINDO/3,  MND0  and  AMI  calculations.  The 
nitrogen  in  the  model  B-lactam  antibiotic  is  predicted  to  invert  periodically  as  the 
dihedral  angle  to  the  exocyclic  N-substituent  sweeps  through  360®. 


0  DISTRIBUTION/ AVAILABILITY  OF  ABSTRACT 
□3  UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT 

□  OTIC  USERS 

21  ABSTRACT  SECURITY  CLASSIFICATION  1 

UNCLASSIFIED  1 

2 a  NAME  OF  RESPONSIBLE  INDIVIDUAL 

JAMES  J.P.  STEWART 

22b  TELEPHONE  (Include  Area  Code) 

(719)472-2655 

22 c  OFFICE  SYMBOL  1 

FJSRL/NC  1 

D  Form  1473,  JUN  86  Previous  editions  are  obsolete  SECURITY  CLASSIFICATION  OF  this  PAGE 

UNCLASSIFIED 


37?  880008 

Numerical  Sensitivity  of  Trajectories  across  Conformational 
Energy  Hypersurfaces  from  Geometry  Optimized  Molecular 
Orbital  Calculations:  AMI,  MNDO,  and  MINDO/3 


Donald  B.  Boyd*  and  David  W.  Smith 

Lilly  Research  Laboratories,  Eli  Lilly  and  Company,  Indianapolis,  Indiana  46285 

James  J.  P.  Stewart 

The  Frank  J.  Seiler  Research  Laboratory,  United  States  Air  Force  Academy,  Colorado  Springs,  Colorado  80840 

Erich  Wimmer 

Cray  Research,  Inc.,  1333  Northland  Drive,  Mendota  Heights,  Minnesota  55120 

Received  6  July  1987;  accepted  22  September  1987 

The  monocyclic  ^-lactam  [[4(S)-methyl-2-oxo-l-azetidinyl]thia]acetic  acid  was  studied  by  the  semi- 
empirical  molecular  orbital  methods  AMI,  MNDO,  and  MINDO/3.  Using  the  reaction  coordinate  option 
in  the  program  MOPAC  on  VAX  and  Cray  X-MP  computers,  the  potential  energy  curve  was  calculated 
for  rotation  of  the  C2-N1-S-C  torsional  angle  in  the  conformationally  flexible  side  chain  while  optimizing 
all  other  geometrical  variables  in  the  molecule.  The  trajectory  taken  during  geometry  optimization  was 
found  to  be  sensitive  to  the  computer,  the  program  version,  the  convergence  criteria,  and  the  degree  of 
code  optimization  used  in  the  calculation.  In  order  to  reduce  the  likelihood  of  spurious  results,  confor¬ 
mational  or  reaction  energy  hypersurfaces  need  to  be  calculated  with  the  more  precise  SCF  convergence 
and  minimization  criteria  available  in  programs  for  MINDO/3,  MNDO,  and  AMI  calculations.  The 
nitrogen  in  the  model  /3-lactam  antibiotic  is  predicted  to  invert  periodically  as  the  dihedral  angle  to  the 
exocyclic  N-substituent  sweeps  through  360°. 


INTRODUCTION 

One  of  the  most  widely  used  computational 
chemistry1-3  programs  currently  available  is 
MOPAC.4, 8  MOPAC  can  run  three  popular 
semiempirical,  all-valence-electron  molecu¬ 
lar  orbital  methods:  MINDO/3,6  MNDO,7  and 
AMI.8  These  methods  are  capable  of  giving 
fairly  reliable  equilibrium  bond  lengths,  bond 
angles,  and  conformations,  molecular  energies 
(heats  of  formation),  ionization  potentials,  di¬ 
pole  moments,  vibrational  frequencies  and 
intensities,  transition  state  geometries,  and 
thermodynamic  properties  of  organic,  phar¬ 
macological,  biochemical,  and  some  inorganic 
molecular  structures.  The  methods  have  been 
reviewed  elsewhere.8-10  MINDO/3,  the  oldest 
of  these  methods,  has  been  used  in  over 
260  papers.11 

The  importance  of  MOPAC  and  kindred 
programs  (MINDO/3,12-14  MNDO,15 16  and 

*To  whom  correspondence  should  be  addressed. 
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AMPAC17)  to  past  and  present  computational 
chemistry  and  computer-assisted  molecular 
design  research  is  exemplified  by  the  fact  that 
a  sizable  proportion  of  all  papers  reporting 
molecular  orbital  results  in  the  chemical 
literature  has  been  based  on  these  programs. 
This  situation  will  continue.  MOPAC  and 
AMPAC  were  the  programs  most  requested 
from  the  Quantum  Chemistry  Program  Ex¬ 
change  (Indiana  University,  Bloomington,  IN) 
in  a  recent  16-month  period.  The  total  num¬ 
ber  of  copies  of  MOPAC  and  AMPAC  dis¬ 
tributed  (230)  far  exceeded  that  for  any  other 
computational  chemistry  program  including 
GAUSSIAN80-UCSF18  (56)  and  MM219,20 
(42).  The  MOPAC  and  AMPAC  programs  are 
sufficiently  comprehensive  and  versatile  that 
they  are  used  both  stand  alone  and  integrated 
with  computer-assisted  molecular  design 
(CAMD)  software  systems.21-23 

MOPAC  is  written4,6  in  transportable 
FORTRAN-77  so  that  it  can  be  run  fairly 
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easily  on  a  number  of  computer  brands,  in¬ 
cluding  the  Digital  Equipment  Corporation 
VAX  computer  series,  which  is  now  almost 
universally  used  in  computer-assisted  mo¬ 
lecular  design  studies.1  However,  as  larger 
and  larger  research  problems  are  attached 
with  MOPAC,  it  is  inevitable  that  super¬ 
computers,  such  as  a  Cray,  will  be  required. 
Already  many  hardware  vendors,  both  major 
and  minor,  have  adapted  or  are  adapting 
MOPAC  to  their  machines. 

We  wish  to  report  our  experience  with 
MOPAC  running  on  a  Cray  X-MP/416.  We 
also  report  central  processing  unit  (CPU) 
times  to  compare  the  efficiency  of  running 
MOPAC  on  Cray  and  VAX  machines. 

An  unexpectedly  large  variability  in  the 
energy  hypersurface  was  found  for  the  test 
case  on  both  the  VAX  and  the  Cray.  This  dis¬ 
covery  suggests  a  vulnerability  of  confor¬ 
mational  predictions  from  this  program  and 
related  programs,  all  of  which  are  based  on 
the  Davidon-Fletcher-Powell  geometry 
optimization  procedure.6, 24-27  The  routine  for 
doing  this  procedure,  FLEPO,  is  used  in  a 
number  of  quantum  mechanical  programs. 
Despite  the  widespread  use  of  MOPAC  and 
AMPAC  and  their  progenitors,  the  minimiza¬ 
tion  technology  used  in  the  programs  has  appar¬ 
ently  been  taken  at  face  value  by  most  users. 

Our  test  case  was  a  pragmatic  one  in¬ 
volving  the  conformation  of  [[4(S)-methyl- 
2-oxo- l-azetidinyl]thia]acetic  acid.  This 
structure  is  a  model  of  the  sulfur  analogue  of 
[[3(S)-(acylamino)-2-oxo-l-azetidinyl]oxy]acetic 
acids  (oxamazins),  a  new  class  of  monocyclic 
P -lactam  antibiotic.28-30  The  reaction  coordi¬ 
nate  was  the  C2-N1-S-C  dihedral  angle  to  the 
S-CH2-COOH  side  chain.  As  recommended  by 
Dewar,10  all  calculations  were  carried  out 
with  full  geometry  optimization,  except,  of 
course,  for  the  dihedral  angle  that  defines  the 
reaction  coordinate.  Thus  a  multidimensional 
potential  energy  surface  (hypersurface)  is  be¬ 
ing  computed.  In  order  to  save  computer  time 
in  generating  trial  density  matrices,  MOPAC 
and  related  programs  use  the  output  geome¬ 
try  from  each  point  in  the  reaction  coordinate 
path  as  input  for  the  next  point.  Thus,  the 
course  of  a  reaction  coordinate  run  can  be 
influenced  by  all  prior  geometries.  This  test 
case  is  obviously  rigorous,  as  well  as  being 
prototypical  of  many  uses  of  MOPAC  in 
computer-assisted  drug  design  studies. 21 29 


H  CH3 
H^-C3— C4— H 
N 

METHODOLOGY 

The  AMI  method,  which  is  the  most  recent 
molecular  orbital  (MO)  method  from  the  Dewar 
group8  and  the  only  method  in  MOPAC  ca¬ 
pable  of  satisfactorily  treating  hydrogen¬ 
bonding  interactions  and  structures  with 
rotational  freedom  about  single  bonds  be¬ 
tween  spz-hybridized  atoms,  was  chosen  for 
most  of  the  calculations.  Optimized  AMI 
parameters  for  sulfur  have  not  yet  been  pub¬ 
lished.  Version  3.1  of  MOPAC  does  AMI  cal¬ 
culations  on  sulfur-containing  structures  us¬ 
ing  sulfur  parameter  values  that  are  identical 
or  similar  to  those  reported  for  the  MNDO 
method.31  While  these  substitute  values  may 
not  be  ideal  for  AMI,  the  MNDO-like  sulfur 
parameters  give  reasonably  good  results  in 
AMI29,32  and  probably  will  not  be  drastically 
different  than  the  ultimate  AMI  values. 
MNDO  and  MINDO/3  are  applied  here  in  ad¬ 
dition  to  AMI.  Regarding  MINDO/3,  it  has 
been  pointed  out33  that  the  sulfur  MINDO/3 
parameters  widely  used  in  MOPAC  and  re¬ 
lated  programs  are  not  those  used  in  the  origi¬ 
nal  study.6  Discrepancies  in  the  MINDO/3 
and  MNDO  parameterization  have  persisted, 
and  the  parameter  values  in  the  programs 
seem  to  have  become  the  de  facto  standard 
through  use.  We  emphasize  that  parameteri¬ 
zation  is  not  at  issue  in  this  investigation.  The 
minimization  procedure  is.  The  vagaries  of 
the  minimization  procedure  could  show  up 
with  any  parameterization. 

The  series  of  conformers  discussed  here  was 
generated  in  a  reaction  coordinate  run  where 
the  dihedral  angle  controlling  the  side  chain 
conformation  is  driven  through  360°  in  30°  in¬ 
crements.  All  the  calculations  used  identical 
geometrical  data  as  input.  All  54  geometrical 
variables,  which  were  specified  in  internal 
coordinates,  had  been  energy  optimized  in  a 
separate  earlier  AMI  calculation  using  the 
default  options  in  MOPAC  Version  3.0  run¬ 
ning  on  a  VAX  8600  computer.29  The  opti¬ 
mized  C2-N1-S-C  dihedral  angle  was  118.5°. 
The  reaction  coordinate  runs  all  started  by 
reoptimizing  the  starting  geometry  with  the 
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C2-N1-S-C  dihedral  angle  at  118.5°  (con- 
former  #1).  The  second  conformer  has  the 
angle  set  at  140°,  the  third  at  170°,  etc.,  until 
at  conformer  #14  the  dihedral  angle  is  again 
at  140°. 

The  floating  point  operations  were  done  in 
double  precision  (64-bit  words)  on  the 
VAX  8600  and  8800  computers.  The  Cray  X- 
MP/416  uses  a  64-bit  single  precision  word, 
but  has  a  shorter  mantissa  (48  bits)  than  that 
of  a  double  precision  word  on  VAX  machines 
(55  bits).  Hence  the  number  of  significant  fig¬ 
ures  in  intermediate-sized  numbers  is  greater 
on  a  VAX,  although  the  Cray  handles  larger 
and  smaller  exponents. 

MOPAC  Versions  3.0  and  3.1  are  the  stan¬ 
dard  ones  available  from  the  Quantum  Chem¬ 
istry  Program  Exchange.4,5  Version  3.11  is  a 
new  version  under  development  by  Stewart34 
for  the  Cray  which  corrects  some  numerical 
instabilities  in  matrix  diagonalization  which 
showed  up  on  that  machine  in  the  quantities 
(1  -q2)m  when  q  approaches  1.  This  arith¬ 
metic,  which  has  been  recoded  in  MOPAC  3.11, 
caused  numerical  errors  of  the  same  order  of 
magnitude  as  the  convergence  criteria.  The 
new  code  also  allows  better  vectorization,  not¬ 
ably  in  the  generation  of  the  Fock  matrix. 

The  complex  convergence  and  optimization 
criteria  built  into  MOPAC  (and  AMPAC)  are 
described  in  the  program  manual  and  docu¬ 
mentation4,5  and  need  not  be  repeated  here. 
In  part  the  default  self-consistent  field  con¬ 
vergence  criterion  is  a  drop  in  electronic  en¬ 
ergy  of  less  than  0.00001  kcal/mol  on  two 
successive  iterations.  For  geometry  optimiza¬ 
tion  the  default  Davidon-Fletcher-Powell 
criteria  can  be  a  projected  change  in  geometry 
of  less  than  0.0001  A  (or  radians)  ("test  on  X”), 
a  projected  decrease  in  energy  of  less  than 
0.001  kcal/mol  ("Herbert’s  test”),  a  gradient 
norm  of  less  than  1.0  kcal  mol-1  A'1  (or 
kcal  mol'1  rad'1) ("test  on  gradient”),  or  a 
change  in  calculated  heats  of  formation  on 
two  successive  iterations  of  less  than 
0.002  kcal/mol  ("heat  of  formation  test”).  By 
using  the  keyword  PRECISE  in  the  input 
data  to  the  program,  the  SCF  criterion  is 
made  100  times  smaller,  and  the  four  geo¬ 
metric  criteria  are  made  100, 100, 500,  and  50 
times  smaller,  respectively.  Depending  on  the 
type  of  calculation  being  done,  the  criteria  are 
adjusted  automatically  by  the  program  to  en¬ 
hance  performance. 


For  the  test  case  studied  here,  all  calcula¬ 
tions  reached  SCF  convergence.  The  number 
of  conformers  which  satisfied  Herbert’s  test 
ranged  from  6  to  14  per  run.  In  the  remaining 
conformers,  the  gradient  test  was  not  passed, 
but  further  work  was  not  justified  because  of 
the  shallowness  of  the  energy  gradient. 

RESULTS 

The  potential  energy  curves  in  Figure  1  are 
labeled  according  to  the  machine  on  which 
they  were  run,  the  version  of  MOPAC  used, 
and  whether  the  calculation  was  done  with 
the  convergence  and  optimization  criteria 
tightened.  Although  all  the  curves  start  out 
similar  and  the  height  of  the  first  barrier 
of  the  twofold  potential  varies  less  than 
0.2  kcal/mol,  some  differentiation  is  seen  at 
the  first  trough  (conformer  #6).  The  height  of 
the  second  barrier  shows  greater  variability 
of  0.3  kcal/mol.  The  energies  of  the  last  three 
conformers  (#12,  #13,  and  #14)  are  highly 
variable,  and,  in  fact,  the  shapes  of  the  curves 
are  obviously  different  qualitatively. 

The  startling  variability  in  the  potential 
energy  curves  in  Figure  1  is  due  to  the  calcu¬ 
lations  following  different  minimum  energy 
pathways  through  the  multidimensional 
space  of  the  geometrical  parameters.  Very 
small  differences  in  the  numerical  accuracy 
of  the  computations  is  causing  the  geometry 
optimization  to  veer  down  diverging  paths 
when  the  potential  energy  surface  offers 
several  channels  of  very  nearly  equal  energy. 
The  effect  is  seen  on  both  the  VAX  and  Cray 
computers. 

The  only  two  curves  of  Figure  1  that  are 
essentially  identical,  as  one  would  hope  that 
hey  all  would  be,  are  for  the  calculations 
done  with  PRECISE.  Although  similar,  the 
heats  of  formation  for  these  two  curves  agree 
only  to  the  first  or  second  decimal  place. 
When  the  calculation  is  done  with  PRECISE 
turned  on,  conformers  #2  and  #14  (both  of 
which  have  dihedral  angle  C2-N1-S-C  =  140°) 
are  more  nearly  identical  in  energy.  Also  with 
PRECISE  the  potential  energy  curves  corre¬ 
spond  more  closely  to  twofold  periodicity. 

Further  tests  were  run  on  the  Cray  X-MP  to 
explore  the  effects  of  using  a  newly  available 
FORTRAN  compiler,  running  the  calculation 
with  single-precision  (64-bit)  or  double¬ 
precision  (128-bit)  words,  and  turning  off  the 
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Figure  1.  AMI  potential  energy  curves  for  [{4(S)-methyl-2-oxo-l-azetidinyl]thia]acetic  acid 
as  predicted  by  MOPAC  on  VAX  and  Cray  computers.  Free  geometry  optimization  was  done 
on  each  conformer  starting  with  an  optimized  geometry.  Each  curve  is  labeled  according  to 
the  machine  and  MOPAC  program  version  (V)  used.  Curves  marked  with  P  were  run  with  the 
PRECISE  option  turned  on  in  MOPAC. 

automatic  vectorization  of  the  code.  The  dif¬ 
ferent  Cray  compilers  achieve  different  levels 
of  automatic  optimization  of  the  operations, 
which  implies  different  orders  of  operations. 

Turning  vectorization  on  or  off  also  implies 
changing  the  order  of  some  operations.  The 
Cray  FORTRAN  compilers  allow  aspects  of 
the  executable  module  to  be  invoked  without 
having  to  change  the  MOPAC  source  code.  As 
seen  in  Figure  2,  these  tests  underscore  the 
sensitivity  of  the  MOPAC  results  when  the 
calculations  are  done  without  PRECISE. 

The  potential  energy  profiles  predicted  by 
MOPAC  using  the  MNDO  and  MINDO/3  mole¬ 
cular  orbital  methods  are  shown  in  Figures  3 
and  4.  The  older  methods  also  show  a  sensi¬ 


tivity  to  the  convergence  and  optimization 
criteria,  although  different  than  AMls.  The 
MINDO/3  curves  start  out  on  a  high-energy 
hypersurface,  but  after  a  few  conformers,  the 
minimizer  finds  a  lower-energy  one.  Subse¬ 
quent  conformers  transverse  the  latter  surface, 
which  has  3-4  kcal/mol  barriers  to  rotation 
about  the  N — S  bond.  The  higher  energy  sur¬ 
face  has  a  gauche  Nl-S-C-C  conformation, 
while  the  lower  one  has  a  trans  conformation 
(see  Table  IV). 

The  "ultimate”  test  was  done  by  running 
the  data  with  PRECISE  in  MOPAC  3.11  in 
double  precision  on  the  Cray  X-MP.  The 
128-bit  words  correspond  to  29  significant 
decimal  places  compared  to  14  for  64-bit 
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Abbreviations: 

P  -  PRECISE  used  0  -  old  1.15  compiler 

V  -  MOPAC  version  N  -  new  1.15  compiler 

S  -  single  precision  0  -  double  precision  - 

M  -  w/o  vectorlzation  Accession  For 
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Figure  2.  AMI  potential  energy  curves  for  ([4(S)-methyl-2-oxo-l-azetidinyl]thia]acetic 
acid  as  predicted  by  MOPAC  on  the  Cray  X-MP/416.  Each  curve  is  labeled  according  to  the 
MOPAC  program  version  (V)  used,  whether  the  calculation  was  done  with  single  (S)  or 
double  (D)  precision,  and  whether  FORTRAN  compilation  was  done  with  the  old  (O)  or  new 
(N)  Releases  1.15  of  the  Cray  CFT  compiler  (see  footnotes  to  Table  I).  The  curve  marked 
W  was  run  without  automatic  vectorization. 


words  normally  used  on  the  Cray.  The  AMI 
method  was  used.  The  curves  in  Figure  5 
show  that  the  calculated  hypersurface  has 
stabilized  for  this  molecule,  and  no  further 
variability  is  revealed  at  the  highest  level  of 
numerical  accuracy. 

The  computer  time  requirements  are  given 
in  Table  I.  MOPAC  ran  the  AMI  calculations 
13-15  times  faster  on  one  processor  of  the 
Cray  X-MP/416  than  on  one  processor  of  a 
VAX  8800.  Vectorization  has  no  effect  on  the 
potential  energy  curve,  and  only  gives  a  two¬ 
fold  improvement  in  speed.  Turning  PRECISE 
on  causes  the  calculations  to  run  2-4  times 
longer.  The  PRECISE,  double-precision  cal¬ 
culation  on  the  Cray  consumed  a  huge  amount 


of  CPU  time,  taking  almost  as  long  as  the 
PRECISE  run  on  the  VAX  8800. 

DISCUSSION  OF  MOLECULAR 
STRUCTURE 

The  seeming  capriciousness  of  the  results 
due  to  small  round-off  errors  being  magnified 
to  chemical  significance  by  the  minimization 
procedure  is  seen  in  the  two  sets  of  Nl-S-C-C 
dihedral  angles  in  Table  II  which  diverge  dra¬ 
matically  at  conformer  #11.  The  differences 
in  the  conformations  can  be  appreciated 
from  the  ORTEP35  ball-and-stick  molecular 
graphics  in  Figure  6. 

The  pyramidal  character  of  the  /3 -lactam 
nitrogen  is  most  clearly  measured  by  the  dis- 


392 


Boyd  et  al.' 


-112-t 


_  -114 
V 


10 

u 

—  -118 


L. 

a 


10 

V 


-iiaw 


-12<H 


-1Z2-I 


-184- 


+  VAX  V3.1 


X  VAX  V3.1  P 


V  -  MOPAC  version 


6 


— r— 
10 


„1 - 1 - T - 1 

11  12  13  14 


Canformer  number 


Figure  3.  MNDO  potential  energy  curves  for  ([4(S)-methyl-2-oxo-l- 
azetidinyl]thia]acetic  acid  as  predicted  by  MOPAC  3.1  on  VAX  8800  and  Cray 
X-MP/416  computers.  The  curves  marked  with  P  were  run  with  the  PRECISE  option 
turned  on. 


tance  h  that  this  atom  is  out  of  the  plane  of  its 
three  substituents.  The  twofold  variation  in 
distance  h  as  a  function  of  the  conformation  of 
the  side  chain  found  with  the  PRECISE  AMI 
calculations  (Table  II)  is  an  interesting  prop¬ 
erty  not  previously  reported.  According  to 
the  minimizer,  the  molecule  reaches  a  lower 
energy  by  inverting  the  /3 -lactam  nitrogen 
when  the  C2-N1-S-C  dihedral  angle  is  near 
- 120°  and  60°.  As  the  side  chain  is  twisted, 
inversion  occurs  upon  approaching  the  mini¬ 
mum  after  the  first  barrier  (Fig.  1).  Then  the 
molecule  inverts  back  to  its  original  config¬ 
uration  when  beginning  to  mount  the  next 
potential  energy  peak.  Inversion  again  occurs 
when  approaching  the  second  minimum,  but 
immediately  flips  back  at  the  next  conformer. 
As  seen  in  Tables  III  and  IV,  the  MNDO  and 


MINDO/3  methods  when  used  with  the  PRE¬ 
CISE  criteria  corroborate  the  inversions  pre¬ 
dicted  by  AMI. 

In  contrast,  in  oxamazin,29  the  oxygen  ana¬ 
logue  of  thiamazin,  the  /3 -lactam  nitrogen  is 
predicted  by  the  AMI  method  (with  PRECISE) 
to  retain  its  configuration  as  the  side-chain 
conformation  changes.  The  oxygen  analogue 
is  predicted  to  have  a  more  pyramidal 
j3 -lactam  nitrogen,29  so  there  is  less  proclivity 
to  invert. 

It  is  also  of  interest  from  a  chemical  point 
of  view  that  AMI  and  MNDO  predict  the 
/3-lactam  nitrogen  of  thiamazin  to  be  much 
more  pyramidal  than  does  MINDO/3 
(Tables  II— IV).  AMI  predicts  the  highest  h 
values  (0.3-0.5  A),  while  MNDO  gives  ex¬ 
tremes  of  0.3-0.4  A.  MINDO/3’s  prediction  of 
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Table  I.  Central  processor  unit  time  required  for  tests  of  MOPAC. 


Machine 

MOPAC 

version 

Cray 

compiler 

PRECISE 

Time  (s) 

Relative 

time 

VAX  8600 

3.0 

AMI 

10888 

28 

VAX  8800 

3.1 

— 

— 

5919 

15 

VAX  8800 

3.1 

— 

on 

20603 

53 

Cray  X-MP/416* 

3.10 

oldb 

— 

393 

1 

Cray  X-MP/416 

3.11 

oldb 

— 

460 

1.2 

Cray  X-MP/416 

3.11 

oidb 

on 

1581 

4 

Cray  X-MP/416 

3.10 

new' 

— 

452 

1.1 

Cray  X-MP/416 

3.10 

new* 

— 

831 

2 

■  Cray  X-MP/416 

3.10 

new' 

— 

7324 

19 

Cray  X-MP/416 

3.11 

new' 

— 

389 

1 

Cray  X-MP/416 

3.11 

newf 

on 

1456 

4 

Cray  X-MP/416 

3.11 

new* 

on 

17950 

46 

VAX  8800 

3.1 

MNDO 

6274 

16 

VAX  8800 

3.1 

— 

on 

16793 

43 

Cray  X-MP/416 

3.11 

new' 

— 

363 

1 

Cray  X-MP/416 

3.11 

new' 

on 

977 

2.5 

VAX  8800 

3.1 

MINDO/3 

3430 

9 

VAX  8800 

3.1 

— 

on 

7595 

20 

Cray  X-MP/416 

3.11 

new' 

— 

157 

0.4 

Cray  X-MP/416 

3.11 

new' 

on 

359 

1 

‘All  calculations  on  the  Cray  were  done  on  one  processor  of  Serial  Number  218,  Mendota  Heights. 

'The  executable  module  was  created  by  the  Cray  compiler  (December  1986-January  1987)  with  default  options  of 
single  precision  and  automatic  vectorization  of  the  code. 

The  executable  module  was  created  by  the  Cray  compiler  CFT  1.15  BF1  (February  1987)  under  COS  1.16  with 
default  options  of  single  precision  and  automatic  vectorization  of  the  code. 

'‘Differs  from  the  default  Cray  executable  by  having  vectorization  turned  off  during  compilation  and  execution. 

‘Differs  from  the  default  Cray  executable  by  having  double-precision  length  (128-bit)  words  for  all  operations. 

fCFT  1.15  BF2  (March  1987)  compiler  under  COS  1.16  using  Cray  libraries  for  RSP  (diagonalization)  routines. 

CCFT  1.15  BF2  (March  1987)  compiler  under  COS  1.16  using  MOPAC’s  internal  RSP  routines  and  with  double- 
precision  length  (128-bit)  words  for  all  operations. 


Table  II.  AMI  optimized  geometrical  data  for  C2-N1-S-C  conformers  of  [[4(S)-methyl-2-oxo-l-azetidinyl]thia]acetic 
acid  obtained  from  MOPAC  3.1  on  a  VAX  8800  without  and  with  PRECISE* 


Conformer 

C2-N1-S-C  (°) 

Nl-S-C-C  O 

h  (A) 

without 

with 

without 

with 

1 

118.5 

-63 

-64 

0.28 

0.27 

2 

140 

-63 

-64 

0.37 

0.38 

3 

170 

-55 

-53 

0.44 

0.45 

4 

-160 

-71 

-73 

0.48 

0.47 

5 

-130 

-77 

-74 

-0.29 

-0.31 

6 

-100 

-76 

-73 

-0.10 

-0.09 

7 

-70 

-76 

-75 

0.14 

0.13 

8 

-40 

-77 

-78 

0.30 

0.30 

9 

-10 

-80 

-82 

0.35 

0.39 

10 

20 

-94 

-94 

0.43 

0.42 

11 

50 

-108 

-117 

0.46 

0.44 

12 

80 

-110 

-52 

0.42 

-0.16 

13 

110 

-109 

-53 

0.39 

0.07 

14 

140 

-109 

-66 

0.41 

0.38 

*A  positive  torsional  angle  for  A-B-C-D  is  measured  clockwise  from  the  A-B-C  plane  to  the  B-C-D  plane  looking 
from  B  to  C.  A  0°  angle  corresponds  to  cis.  h  is  the  distance  the  /3-lactam  nitrogen  is  out  of  the  plane  of  its  three 
substituents.  A  negative  value  of  h  indicates  sulfur  is  on  the  same  side  of  the  amide  plane  as  the  4a-methyl  carbon. 
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Figure  4.  MINDO/3  potential  energy  curves  for  [[4(S)-methyl-2-oxo-l- 
azetidinyl]thia]acetic  acid  as  predicted  by  MOPAC  3.1  on  VAX  8800  and  Cray 
X-MP/416  computers.  The  curves  marked  with  P  were  run  with  the  PRECISE 
option  turned  on  in  MOPAC. 


Table  III.  MNDO  optimized  geometrical  data  for  C2-N1-S-C  conformers  of  [[4(S)-methyl-2-oxo-l 
azetidinyl}thia]acetic  acid  obtained  from  MOPAC  3.1  on  a  VAX  8800  without  and  with  PRECISE. 


Conformer 

C2-N1-S-C  n 

Nl-S-C-C  (°) 

h  (A) 

without 

with 

without 

with 

1 

118.5 

-64 

-66 

0.28 

0.29 

2 

140 

-64 

-65 

0.36 

0.36 

3 

170 

-fr 

-68 

0.43 

0.43 

4 

-160 

-83 

-89 

0.40 

0.38 

5 

-130 

-89 

-95 

-0.27 

-0.28 

6 

-100 

-88 

-92 

-0.15 

-0.13 

7 

-70 

-87 

-88 

-0.03 

0.04 

8 

-40 

-85 

-79 

0.18 

0.18 

9 

-10 

-85 

-76 

0.24 

0.25 

10 

20 

-94 

-81 

0.26 

0.25 

11 

50 

-110 

-88 

0.14 

-0.07 

12 

80 

-111 

-87 

0.07 

0.08 

13 

110 

-109 

-85 

0.25 

0.26 

14 

140 

-105 

-83 

0.37 

0.37 
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Figure  5.  AMI  potential  energy  curves  for  [[4(S)-methyl-2-oxo-l- 
azetidinyl  lthia)acetic  acid  as  predicted  by  MOPAC  on  VAX  8800  and  Cray  X-MP/416 
computers.  All  were  computed  wivh  the  PRECISE  option  turned  on.  Curves  are 
labeled  as  in  Figure  2. 


Table  IV.  MINDO/3  optimized  geometrical  data  for  C2-N1-S-C  conformers  of  [[4(S)-methyl-2  oxo-1- 
azetidinyl]thia)acetic  acid  obtained  from  MOPAC  3.1  on  a  VAX  8800  without  and  with  PRECISE. 


Conformer 

C2-N1-S-C  (°) 

Nl-S-C-C  (°) 

h  (A) 

without 

with 

without 

with 

1 

118.5 

-65 

-68 

0.14 

0.14 

2 

140 

-67 

-79 

0.18 

0.17 

3 

170 

-88 

180 

0.14 

0.12 

4 

-160 

179 

180 

-0.09 

-0.09 

5 

-130 

179 

179 

-0.10 

-0.10 

6 

-100 

179 

178 

-0.05 

-0.04 

7 

-70 

179 

178 

0.02 

0.02 

8 

-40 

179 

180 

0.05 

0.06 

9 

-10 

180 

180 

0.04 

0.05 

10 

20 

180 

180 

0.03 

0.01 

11 

50 

180 

180 

-0.01 

-0.02 

12 

80 

180 

-179 

0.02 

0.04 

13 

110 

-179 

-178 

0.12 

0.12 

14 

140 

-179 

-  178 

0.17 

0.17 
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Figure  6.  Conformers  #2  (top),  #12  (middle),  and  #14 
(bottom)  from  the  AMI  calculation  done  using  MOPAC 
Version  3.1  on  the  VAX  8800  without  (left)  and  with 
(right)  PRECISE.  Conformers  #2  and  #14  both  have  a 
C2-N1-S-C  dihedral  angle  of  140°,  and  the  sulfur  and 
4a  -methyl  carbon  are  on  opposite  sides  of  the  amide 
plane.  The  /J-lactam  nitrogen  has  inverted  in  conformer 
#12  depending  on  whether  PRECISE  is  used  in  the 
MOPAC  reaction  coordinate  nm.  The  0 -lactam  car¬ 
bonyl  oxygen  projects  toward  the  viewer. 


0.14  A  comes  closest  to  the  experimental  h 
value  of  0.18  A  for  a  thiamazin  with  a 
C2-N1-S-C  torsional  angle  at  119°  in  the  solid 
state.29 

The  distance  h  has  been  associated  with 
antibacterial  activity  of  /3 -lactam  compounds.36 
Following  Sweet’s  crystallographic  structure- 
activity  relationship  studies,  it  became  com¬ 
monly  believed  that  compounds  with  higher  h 
values  have  better  biological  activity.  Accu¬ 
mulating  crystallographic  data37,38  have 
shown  that  this  concept  is  not  as  universally 
true  as  once  thought. 

From  the  point  of  a  selecting  a  computa¬ 
tional  chemistry  method  that  will  give  re¬ 
liable  h  values,  the  venerable  MINDO/3 
method  is  actually  superior  to  the  more  pol¬ 
ished  and  generally  useful  MNDO  and  AMI 
methods.  Also  MINDO/3  predicts  the  equi¬ 
librium  -lactam  C  —  N  bond  length  to  be 
1.37-1.38  A  (depending  on  side-chain  con¬ 
formation),  which  compares  favorably  with 
the  observed  X-ray  diffraction  value29  of 
1.39  A.  In  contrast,  AMI  overestimates  this 
length  at  1.42-1.44  A,  and  MNDO  over¬ 
estimates  it  even  more  at  1.43-1.45  A. 


Corresponding  to  the  twofold  potentials 
(Fig.  1-5),  the  calculated  preferred  confor¬ 
mations  of  the  model  structure  are  gauche 
and  thus  compatible  with  the  crystalline 
state  conformation  of  a  thiamazin.29  The  two 
low-energy  conformers  are  near  -100°  and 
110°.  Indications  from  temperature-dependent" 
nuclear  magnetic  resonance  (NMR)  studies 
are  that  sulfenamides  in  general  have  bar¬ 
riers  to  rotation  around  the  N — S  bond  of 
9-23  kcal/mol.39  Acylation  of  the  nitrogen  of 
a  /3-lactam  would  lower  this  due  to  delocaliza¬ 
tion  of  the  nitrogen  lone  pair.  Nevertheless 
the  barrier  heights  of  4-6  kcal/mol  (Fig.  1) 
that  are  obtained  with  the  unoptimized  AMI 
sulfur  parameters  appear  low.29 

Even  though  all  three  of  the  molecular 
orbital  methods6"8  were  parameterized  by 
Dewar  and  coworkers  on  the  basis  of  repro¬ 
ducing  experimental  heats  of  formation  for  a 
wide  variety  of  molecules,  the  variety  did  not 
include  azetidinones.  Consequently,  it  is  not 
surprising  that  the  predicted  AHfs  for  our 
/3 -lactam  structure  vary  over  a  wide  range 
(Fig.  1, 3,  and  4).  AMI  gives  AHf  s  ca.  -80  kcal/ 
mol,  MNDO  is  intermediate  at  ca.  -120  kal/ 
mol,  and  MINDO/3,  which  is  known  to  be 
poor  for  four-membered  rings,3  gives  ca. 
- 180  kcal/mol. 

CONCLUSION 

The  conclusion  from  these  tests  is  that 
MOPAC  when  run  with  the  current  default 
options  is  not  as  reliable  in  terms  of  repro¬ 
ducibility  as  many  users  may  have  assumed. 
When  certain  conformational  energy  hyper¬ 
surfaces  are  being  traversed,  the  molecule 
can  go  along  different  paths  on  the  surface 
depending  on  the  computer,  version,  program 
options,  etc.  The  numerical  accuracy  of  the 
calculations  is  having  a  significant  effect  on 
the  apparent  shape  of  the  surface. 

There  is  no  way  of  estimating  how  many 
published  calculations  of  transition  state 
geometries  and  conformations  may  have  been 
affected  by  this  hypersensitivity  of  molecular 
orbital  programs  based  on  the  Davidon- 
Fletcher-Powell24-27  geometry  optimization 
procedure.  Very  few  papers  in  the  literature 
reporting  MINDO/3,  MNDO,  and  AMI  re¬ 
sults  have  been  explicit  about  whether  the 
calculations  were  done  with  the  tighter  SCF 
convergence  and  minimization  criteria.  We 
call  upon  authors  to  be  more  specific  in 
the  future. 
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The  problem  with  the  minimization  proce¬ 
dure  uncovered  here  will  not  occur  for  every 
molecule.  It  is  most  likely  to  occur  when  the 
hypersurface  traversed  is  flat.  However,  it 
would  be  difficult  to  know  ahead  of  time  when 
such  a  flat  surface  will  be  encountered  for  a 
complex  molecule.  Hence  the  only  safe  course 
is  to  always  ny  out  and/or  use  the  tightest 
available  convergence  and  minimization 
criteria.  As  seen  in  Figure  5,  use  of  PRECISE 
helps  guarantee  that  the  potential  energy 
curves  are  at  least  visually  similar. 

In  order  to  make  MOPAC  as  reliable  as 
possible,  a  new  version  with  a  superior  minimi¬ 
zation  method  is  under  development  and  will 
be  released  after  testing.41  It  will  also  have 
enhancements  stemming  from  experience  of 
running  MOPAC  on  the  Cray.  Until  the  new 
version  is  available,  calculations  with  MOPAC, 
AMPAC,  and  consanguineous  programs 
should  be  done  with  the  tightest  available 
convergence  anu  minimization  criteria.  Such 
a  procedure  will  cost  several  times  more  com¬ 
puter  time  (see  Table  I),  but  such  costs  are 
necessary  to  help  insure  stable  results. 

For  the  test  case40  discussed  in  this  paper, 
MOPAC  ran  13-21  times  faster  on  the  Cray 
X-MP/416  than  on  a  VAX  8800.  This  in¬ 
crease  in  speed  is  encouraging,  and  even 
greater  increases  in  speed  may  be  possible  if 
the  code  were  fully  customized  for  the  Cray. 
The  speed  enhancements  achievable  with  a 
supercomputer  will  permit  more  ambitious 
computational  experiments  to  be  performed 
even  with  fine  convergence  and  optimization 
criteria. 

One  of  us  (JJPS)  thanks  the  National  Research  Coun¬ 
cil  for  an  award  of  an  Associateship.  Mr.  John  Lacher 
assisted  with  telecommunications,  Ms.  Kathleen  A. 
Cloutier  assisted  with  literature  retrieval,  and  Mr. 
Richard  W.  Counts  furnished  data  about  the  quantities 
of  programs  distributed  by  QCPE. 
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